# -*- coding: utf-8 -*

'''
Created on 23/11/2013

@author: nicolas
'''

import math
import pylab as p  
from scipy import integrate
from bpsk_modular_bits import modular_bits_en_bpsk, generar_senal  


def obtener_muestras_de_entrenamiento(N):
    # genera vector aleatorio de enteros entre 0 y 1
    B = p.random_integers(0, 1, N)
    return B.tolist()



def generar_bit_Bernoulli(pb):
    # tomo un valor aleatorio entre 0 y 1, si es <= que p retorna 1, 0 en caso contrario
    rand = p.rand()
    if (rand <= pb):
        return 1
    else:
        return 0

def generar_bits(lambda_01, lambda_11, N):
    B = []
    # Agrego el primer elemento como una Bernoulli(lambda_01)
    B.append(generar_bit_Bernoulli(lambda_01))
    
    n = 1
    # Siguiendo el esquema de la cadena de markov, agrego los siguientes bits
    while (n < N):
        n = n + 1
        # B[-1] = ultimo elemento de la lista
        if(B[-1] == 0):
            B.append(generar_bit_Bernoulli(lambda_01))
        else:
            B.append(generar_bit_Bernoulli(lambda_11))

    return B



def sumar_ruido_gaussiano_blanco(X, sigma2_W):
    # genero vector de ruido gaussianno
    ruido = p.normal(0, sigma2_W, len(X))
    
    # se lo sumo a la senal X
    X_R = p.sum([X, ruido], 0)      
      
    return X_R

def sumar_ruido_gaussiano_correlacionado(X, sigma2_W):
    # genero vector de ruido gaussianno
    ruido = p.array([0.0]*len(X))
    
    for i in range(len(X)):
        ruido[i] = p.normal(-X[i]*0.25, sigma2_W, 1)
        
    # se lo sumo a la senal X
    X_R = p.sum([X, ruido], 0)      
      
    return X_R

#recive un angulo y devuelve el equivalente entre -pi/2 y 3/2pi
def normalizar_alpha(alpha):
    while (alpha < -1 *math.pi /2):
        alpha = alpha + 2 * math.pi
    while (alpha >  3* math.pi / 2):
        alpha = alpha - 2 * math.pi
    return alpha
        

def estimar_fase_min(S_R, L):
    #encuentro el n que hace que S_R[n] sea minimo
    n_min = p.argmin(S_R)
    alpha_i = (n_min + 1) * 2 * math.pi / L
    alpha_i = normalizar_alpha(alpha_i)
    return [alpha_i, n_min]

def estimar_fase_max(S_R, L):
    #encuentro el n que hace que S_R[n] sea maximo
    n_max = p.argmax(S_R)
    alpha_i = (n_max - L / 2.0 + 1) * 2 * math.pi / L
    alpha_i = normalizar_alpha(alpha_i)
    return [alpha_i, n_max]

def obtener_alpha_mas_proximo_amplitud(S_R, alpha1, alpha2, A):
    err1 = A**2 - p.absolute(S_R[alpha1[1]]) 
    err2 = A**2 - p.absolute(S_R[alpha2[1]])
    if (p.absolute(err1) < p.absolute(err2)):
        return alpha1[0]
    else:
        return alpha2[0]


def estimar_fase_de_senial_bpsk(S_R, A):
    senal = generar_senal(A, len(S_R))
    
    Y = p.convolve(S_R, senal)
    
    alpha_i_max = estimar_fase_max(Y,len(S_R))
    alpha_i_min = estimar_fase_min(Y,len(S_R))
    
    # estimo fase por maximo y minimo
    # si son iguales lo retorno, si no, retorno el que mas se aproxime en modulo a la amplitud
    if (p.absolute(alpha_i_max[0] - alpha_i_min[0]) < 0.01):
        return alpha_i_max[0]
    else:
        return obtener_alpha_mas_proximo_amplitud(Y, alpha_i_max, alpha_i_min, A)

def encontrar_alphas(entrenamiento, fases):
    array_alpha_0 = []
    array_alpha_1 = []
    
    #para cada elemento de las fases veo si corresponde a un 0 o 1 y lo separo en dos vectores
    for i in range(len(entrenamiento)):
        if (entrenamiento[i] == 0):
            array_alpha_0.append(fases[i])
        else:
            array_alpha_1.append(fases[i])
            
    #los valores a usar son el promedio de los estimados
    alpha_0 = p.average(array_alpha_0)
    alpha_1 = p.average(array_alpha_1)
    
    return [alpha_0, alpha_1]

def Q(z):
    result = integrate.quad(lambda x: (1 / math.sqrt((2 * math.pi))) * (math.e ** ((-1) * (x ** 2) / 2)), z, p.Inf)
    return result[0]

def detectar_bayesianamente(S_R, A, alpha_0, alpha_1, C_W, P_0, P_1, P_01_function=None):
    C_W_inv = p.inv(C_W)
    S_R_trp = p.transpose(S_R)
    
    S_0 = p.array(modular_bits_en_bpsk([0], alpha_0, alpha_1, A, len(S_R)))
    S_1 = p.array(modular_bits_en_bpsk([1], alpha_0, alpha_1, A, len(S_R)))
    S_10 = S_1 - S_0

    S_R0 = ((S_0 + S_1) / 2) + (p.log(P_0 / P_1) * S_10) / (p.dot(p.dot(p.transpose(S_10) , C_W_inv) , S_10))
    N_trp = p.transpose(p.dot(C_W_inv , (S_0 - S_1)))

    # regla
    test = (p.dot(p.dot(S_R_trp, C_W_inv) , S_1)) - (p.dot(p.dot(S_R_trp, C_W_inv) , S_0))
    
    B = 0 if  test < p.dot(N_trp , S_R0) else 1
       
    # probabilidad de error     
    r = p.sqrt(p.dot(p.dot(p.transpose(S_10) , C_W_inv) , S_10))
    P_e = (Q((p.log(P_0 / P_1) / r) + (r / 2)) * P_0) + (Q((-p.log(P_0 / P_1) / r) + (r / 2)) * P_1)
    
    return [B, P_e]

def calcular_estado_markov(lambda_01, lambda_11, estado_inicial, n):
    M_markov = p.mat([[1 - lambda_01, lambda_01], [1 - lambda_11, lambda_11]])
    
    # TODO: aca no se pq me esta devolviendo un vector columna, probe de trasponerlo pero no logre arreglarlo.
    # lo transformo en vector fila
    estado_tmp = p.dot(estado_inicial, p.matrix_power(M_markov, n))
    estado = p.array([estado_tmp[0, 0], estado_tmp[0, 1]])
    return estado

def calcular_estado_markov_estacionario(lambda_01, lambda_11):
    p1 = (1 - lambda_11) / (lambda_01 + (1 - lambda_11))
    return p.array([p1, 1 - p1])
